setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210120.H3k.analysis/117012.analysis")
library(data.table)
part1=fread("part1/part1_zscore.tsv",head=T,sep="\t")
names(part1)=c("chrom",names(part1)[2:2010])
part1$unitID=paste(part1$chrom,part1$pos,part1$ref,part1$alt,sep=":")

part2=fread("part2/part2_zscore.tsv",head=T,sep="\t")
names(part2)=c("chrom",names(part2)[2:2010])
part2$unitID=paste(part2$chrom,part2$pos,part2$ref,part2$alt,sep=":")

part3=fread("part3/part3_zscore.tsv",head=T,sep="\t")
names(part3)=c("chrom",names(part3)[2:2010])
part3$unitID=paste(part3$chrom,part3$pos,part3$ref,part3$alt,sep=":")

part4=fread("part4/part4_zscore.tsv",head=T,sep="\t")
names(part4)=c("chrom",names(part4)[2:2010])
part4$unitID=paste(part4$chrom,part4$pos,part4$ref,part4$alt,sep=":")

part5=fread("part5/part5_zscore.tsv",head=T,sep="\t")
names(part5)=c("chrom",names(part5)[2:2010])
part5$unitID=paste(part5$chrom,part5$pos,part5$ref,part5$alt,sep=":")

part6=fread("part6/part6_zscore.tsv",head=T,sep="\t")
names(part6)=c("chrom",names(part6)[2:2010])
part6$unitID=paste(part6$chrom,part6$pos,part6$ref,part6$alt,sep=":")

deepsea =rbind(part1,part2,part3,part4,part5,part6)
deepsea =as.data.frame(deepsea)
rm(part1)
rm(part2)
rm(part3)
rm(part4)
rm(part5)
rm(part6)

###进行一致性的分析
coln=names(deepsea)
braincell=c("Brain_Angular_Gyrus","Brain_Anterior_Caudate","Brain_Cingulate_Gyrus","Brain_Germinal_Matrix","Brain_Hippocampus_Middle","Brain_Inferior_Temporal_Lobe","Brain_Mid_Frontal_Lobe","Brain_Substantia_Nigra")

file1=read.csv("./117K.ASH.alt.up.or.down.statis.csv",header=T)
file2=read.csv("./117K.ASH.add.enh.promtr.csv",header=T)
names(file2)=c("unitID","snp.location")
names(file1)=c("unitID",names(file1)[2:5])

file=merge(file1,file2,by="unitID")
locations=as.character(unique(file$snp.location))

for(lc in locations){
ASH =file[file$snp.location==lc,]#按位置进行分类
for(i in 1:length(braincell)){
braincn=coln[grep(coln,pattern = braincell[i])]
braincn=braincn[grep(braincn,pattern = "H3K")]#筛出H3k的列名
deepseq= deepsea[,c("unitID",braincn)]
result=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
H3k=unlist(strsplit(braincn,"\\|"))[seq(2,3*length(braincn),3)]
colnames(result)=H3k
frt=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
colnames(frt)=H3k
deepseq=merge(deepseq,ASH,by="unitID")#经查验，merge后的deepseq的位点排列方式与ASH一致
for(j in 1:length(braincn)){
result[,j]=ifelse(deepseq[,1+j]<0,"down","up")
frt[,j]=ifelse(result[,j]==deepseq$alt.group,"same","opposite")
}
rt_statis=data.frame(matrix(NA,length(braincn),4))
row.names(rt_statis)=H3k
colnames(rt_statis)=c("same","opposite","same_ratio","pvalue")
for(k in 1:length(braincn)){
rt_statis[k,1]=table(frt[,k]=="same")[2]
rt_statis[k,2]=table(frt[,k]=="opposite")[2]
rt_statis[k,3]=rt_statis[k,1]/(rt_statis[k,1]+rt_statis[k,2])
rt_statis[k,4]=binom.test(rt_statis[k,1],rt_statis[k,1]+rt_statis[k,2],p=0.5)$p.value
}
fn=paste0("./result/consider.all.beta0.",braincell[i],".",lc,"_statis.csv")
write.csv(rt_statis,fn,quote=F,row.names = T)
}
	#算均值然后求差异
sel1=paste0("./result/consider.all.beta0.",braincell,".",lc,"_statis.csv")
i=1
rt=read.csv(sel1[i],header = T)[,1:3]
for (i in 2:length(braincell)) {
  f1=read.csv(sel1[i],header = T)[,1:3]
  rt=rbind(rt,f1)
}

H3k.group=unique(rt$X)
col_names=c("H3k.group","same.mean","opposite.mean","same.ratio","P.value")
result=data.frame(matrix(NA,1,ncol=5))
names(result)=col_names
result=result[-1,]
for (i in 1:length(H3k.group)) {
  resultmp=data.frame(matrix(NA,1,ncol=5))
  names(resultmp)=col_names
  tmp=rt[rt$X==H3k.group[i],]
  resultmp[1,1]=tmp[1,1]
  resultmp[1,2]=round(mean(tmp[,2]),digits = 0)
  resultmp[1,3]=round(mean(tmp[,3]),digits = 0)
  resultmp[1,4]=resultmp[1,2]/(resultmp[1,2]+resultmp[1,3])
  resultmp[1,5]=binom.test(resultmp[1,2],(resultmp[1,2]+resultmp[1,3]),p=0.5)$p.value
  result=rbind(result,resultmp)
}
fn2=paste0("./result/117K.ASH.",lc,".statis.csv")
write.csv(result,fn2,quote=F,row.names = F)
}


#ASM部分

asm=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/869727.all.snp.vaf.up.down",head=T,sep="\t")
asm2=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/220520ASMs_anno.hg19_multianno.csv",head=T,sep=",")	#220K
asm2$unitID=paste(asm2$Chr,asm2$Start,asm2$Ref,asm2$Alt,sep=":")
asm=asm[asm$rt.unitID %in% asm2$unitID,]
filea=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
asm=asm[asm$rt.unitID %in% filea$unitID,]
asm=merge(file2,asm,by.y = "rt.unitID",by.x="unitID")

dim(asm)
[1] 13649     6

locations =as.character(unique(asm$snp.location))
coln=names(deepsea)
for(lc in locations){
ASH =asm[asm$snp.location==lc,]#按位置进行分类
for(i in 1:length(braincell)){
braincn= coln[grep(coln,pattern = braincell[i])]
braincn=braincn[grep(braincn,pattern = "H3K")]
deepseq= deepsea[,c("unitID",braincn)]
result=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
H3k=unlist(strsplit(braincn,"\\|"))[seq(2,3*length(braincn),3)]
colnames(result)=H3k
frt =data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
colnames(frt)=H3k
deepseq=merge(deepseq,ASH,by="unitID")
for(j in 1:length(braincn)){
result[,j]=ifelse(deepseq[,1+j]<0,"down","up")
frt[,j]=ifelse(result[,j]==deepseq$group.up.down,"same","opposite")
}
rt_statis=data.frame(matrix(NA,length(braincn),4))
row.names(rt_statis)=H3k
colnames(rt_statis)=c("same","opposite","same_ratio","pvalue")
for(k in 1:length(braincn)){
rt_statis[k,1]=table(frt[,k]=="same")[2]
rt_statis[k,2]=table(frt[,k]=="opposite")[2]
rt_statis[k,3]=rt_statis[k,1]/(rt_statis[k,1]+rt_statis[k,2])
rt_statis[k,4]=binom.test(rt_statis[k,1],rt_statis[k,1]+rt_statis[k,2],p=0.5)$p.value
}
fn=paste0("./result.ASM.220K.ovarlap.117K/consider.all.beta0.",braincell[i],".",lc,"_statis.csv")
write.csv(rt_statis,fn,quote=F,row.names = T)
}
	#算均值然后求差异
sel1=paste0("./result.ASM.220K.ovarlap.117K/consider.all.beta0.",braincell,".",lc,"_statis.csv")
i=1
rt=read.csv(sel1[i],header = T)[,1:3]
for (i in 2:length(braincell)) {
  f1=read.csv(sel1[i],header = T)[,1:3]
  rt=rbind(rt,f1)
}

H3k.group=unique(rt$X)
col_names=c("H3k.group","same.mean","opposite.mean","same.ratio","P.value")
result=data.frame(matrix(NA,1,ncol=5))
names(result)=col_names
result=result[-1,]
for (i in 1:length(H3k.group)) {
  resultmp=data.frame(matrix(NA,1,ncol=5))
  names(resultmp)=col_names
  tmp=rt[rt$X==H3k.group[i],]
  resultmp[1,1]=tmp[1,1]
  resultmp[1,2]=round(mean(tmp[,2]),digits = 0)
  resultmp[1,3]=round(mean(tmp[,3]),digits = 0)
  resultmp[1,4]=resultmp[1,2]/(resultmp[1,2]+resultmp[1,3])
  resultmp[1,5]=binom.test(resultmp[1,2],(resultmp[1,2]+resultmp[1,3]),p=0.5)$p.value
  result=rbind(result,resultmp)
}
fn2=paste0("./result.ASM.220K.ovarlap.117K/117K.ASH.overlap.ASM.",lc,".statis.csv")
write.csv(result,fn2,quote=F,row.names = F)
}
